#include "navierstokes.hpp"

// variational formulation language
#include <life/lifevf/vf.hpp>

namespace Life
{


template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void
NavierStokes<Dim, Order, Entity>::updateResidual2( const vector_ptrtype& X, vector_ptrtype& R )
{

    using namespace Life::vf;

    mesh_ptrtype mesh = Xh->mesh();

    element_fluid_type U( Xh, "u" ); U = *X;
    element_fluid_type V( Xh, "v" );
    element_fluid_0_type u = U.template element<0>();
    element_fluid_0_type v = V.template element<0>();
    element_fluid_1_type p = U.template element<1>();
    element_fluid_1_type q = V.template element<1>();


    AUTO( defv, 0.5*(gradv(u)+trans(gradv(u)) ) );
    AUTO( def, 0.5*(grad(v)+trans(grad(v)) ) );

    //Matrice Identitée
    AUTO (Id , oneX()*trans(oneX()) + oneY()*trans(oneY()) );

    // Tenseur des contraintes (trial)
    AUTO(Sigmav, -idv(p)*Id + 2*mu*defv );

    // Tenseur des contraintes (test)
    AUTO(Sigma, -id(q)*Id + 2*mu*def );

    //force volumique : 
    AUTO(f,vec(cst(0.),cst(0.)/*-rho*cst(9.81)*/ ) );

    // Vitesse d'entrée du fluide
    double vvv;
    if ( M_bdf_fluid->time()<=0.5) { vvv=M_bdf_fluid->time()/0.5; }
    else {vvv=1.0;}

    int stationnaire = this->application()->vm()["stat"].template as<int>();
    if (stationnaire==1) vvv=Vitesse_Entree;

    AUTO(Vit_Entree, vvv*Vitesse_Entree*( Py()*(1.0-Py()) )*oneX() );


   //--------------------------------------------------------------------------------------------------//
#if 0
    //le terme non linéaire
    form1( Xh, R, _init=true ) =
        integrate( elements(mesh), _Q<3*(Order)>(),
              - rho*trans(idv(u))*gradv(u)*id(v) );
#endif
    // sigma : grad(v) sur Omega
    form1( Xh, R ) += 
        integrate( elements(mesh), _Q<2*(Order)>(),
            - trace(Sigmav*trans(grad(v))) ); 
   
    // terme assurant l'incompressibilité
    form1( Xh, R ) += 
        integrate( elements(mesh), _Q<2*(Order-1)>(),
            - divv(u)*id(q) );

   //--------------------------------------------------------------------------------------------------//
#if 1
    // terme de bord classique
    form1( Xh, R ) +=
        integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
            + trans(Sigmav*N())*id(v) );
    
    // terme de bord classique    
    form1( Xh, R ) +=
        integrate( markedfaces(mesh,mesh->markerName("Paroi")), _Q<2*Order>(),
            + trans(Sigmav*N())*id(v) );

    //terme de bord weak boundary
    form1( Xh, R ) +=
        integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
            + trans(Sigma*N())*idv(u)
            - penalbc*trans(idv(u))*id(v)/hFace() );

    //terme de bord weak boundary
    form1( Xh, R ) +=
        integrate( markedfaces(mesh,mesh->markerName("Paroi")), _Q<2*Order>(),
            + trans(Sigma*N())*idv(u)
            - penalbc*trans(idv(u))*id(v)/hFace() );     
#endif
   //--------------------------------------------------------------------------------------------------//

    //Force Volumique
    form1( Xh, R ) +=
        integrate( elements(mesh), _Q<2*Order>(),
            + trans(f)*id(v) );

   //--------------------------------------------------------------------------------------------------//

#if 1
    //terme de bord weak boundary
    form1( Xh, R ) +=
        integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<3*Order>(),
            - trans(Sigma*N())*Vit_Entree //+ trans(2*mu*def*N())*Vit_Entree//- trans(Sigma*N())*Vit_Entree
            + penalbc*trans(Vit_Entree)*id(v)/hFace() );
#endif

/*
 // terme de bord classique    
    form1( Xh, R ) +=
        integrate( markedfaces(mesh,mesh->markerName("Sortie")), _Q<2*Order>(),
            //trans(+idt(p)*N()-mu*deft*N())*id(v)); //- trans(Sigmat*N())*id(v) ); ///////////////////
            - trans(Sigmav*N())*id(v) );

    form1( Xh, R ) +=
        integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
            + penalbc*idv(p)*id(q)/hFace() );

    form1( Xh, R ) +=
        integrate( markedfaces(mesh,mesh->markerName("Sortie")), _Q<2*Order>(),
            + penalbc*idv(p)*id(q)/hFace() );

    form1( Xh, R ) +=
        integrate( markedfaces(mesh,mesh->markerName("Sortie")), _Q<2*Order>(),
            - penalbc*cst(100.0)*id(q)/hFace() );
*/
   //--------------------------------------------------------------------------------------------------//

if (stationnaire==0) {
    form1( Xh, R ) +=
        integrate( elements(mesh), _Q< 2*Order >(),
            - trans(idv(u))*id(v)*M_bdf_fluid->polyDerivCoefficient(0)
            + trans(idv( M_bdf_fluid->polyDeriv().template element<0>() ) ) *id(v) );
}

 int stabilisation = this->application()->vm()["stab"].template as<int>();
if (stabilisation==1) {
    //terme de stabilisation
    form1( Xh, R ) +=
        integrate( internalfaces(mesh), _Q< 3*Order >(),
           // + 0.025*vf::pow(hFace(),2.0)*( abs( trans(idv(u))*N() ))*(trans(jumpv(gradv(u)))*jump(grad(v)) )  );
            + (0.05+vf::pow(hFace(),2.0))*( (trans(idv(u))*N())*( trans(idv(u))*N() ))*(trans(jumpv(gradv(u)))*jump(grad(v)) )  );        
}
   
 
    R->close();

}

} // Life


// instantiation
template class Life::NavierStokes<2,2,Life::Simplex>;
//template class Life::NavierStokes<3,2,Life::Simplex>;
